1.3.1 Software

• Installieren Sie R und seine Freunde.
• Installieren Sie die folgende R-Pakete:
  • tidyverse
  • tidymodels
  • weitere Pakete werden im Unterricht bekannt gegeben (es schadet aber nichts, jetzt schon Pakete nach eigenem Ermessen zu installieren)
• R Syntax aus dem Unterricht findet sich im Github-Repo bzw. Ordner zum jeweiligen Semester.

1.3.2 Videos

• Playlist zu den Themen
• Auf dem YouTube-Kanal des Autors finden sich eine Reihe von Videos mit Bezug zum Inhalt dieses Buches.

1.3.3 Online-Zusammenarbeit

Hier finden Sie einige Werkzeuge, die das Online-Zusammenarbeiten vereinfachen:

• Frag-Jetzt-Raum zum anonymen Fragen stellen während des Unterrichts. Der Keycode wird Ihnen vom Dozenten bereitgestellt.
• Padlet zum einfachen (und anonymen) Hochladen von Arbeitsergebnissen der Studentis im Unterricht. Wir nutzen es als eine Art Pinwand zum Sammeln von Arbeitsbeiträgen. Die Zugangsdaten stellt Ihnen der Dozent bereit.

2.7.1 Abschnitt Forschungsfrage und Hintergrund

In diesem Abschnitt passiert noch keine Statistik bzw. keine Analyse. Stattdessen stellen Sie in “normaler Sprache”, also ohne intensiven Gebrauch vom (statistischem) Fachvokabular dar, was Ziel und was Hintergrund der Analyse ist. Sie können als Ziel bzw. Hintergrund den formalen Aspekt der Prüfung anführen, wichiger sind aber inhaltliche bzw. fachliche Überlegungen: Worum geht es in der Analyse? Warum ist die Frage wichtig? Was wird untersucht? Anhand welcher Methodik wird die Frage untersucht?

2.7.2 Vorbereitung

In diesem Abschnitt Ihres Analysedokuments führen Sie die technische Vorbereitung durch. Das betrifft vor allem das Importieren der Daten und das Starten aller R-Pakete, die in der Analyse verwendet werden.

2.7.3 Explorative Datenanalyse

Die explorative Datenanalyse (EDA) meint sowohl die deskriptive Statistik als auch die Datenvisualisierung. Typische Schritte sind: das Bearbeiten (oder Entfernen) von Extremwerten und fehlenden Werten, die Untersuchung von Verteilungsformen oder das Suchen nach Mustern (Korrelationen, Gruppenunterschieden). Ein nützliches Ergebnis ist z.B. zu erkennen, welche Variablen sich als Prädiktoren eignen (für den nächsten Abschnitt der Modellierung). Ziel ist, dass Sie den folgenden Schritt vorbereiten, also Schritte unternehmen, damit Sie die AV möglichst gut vorhersagen können.

2.7.4 Modellierung

In diesem Schritt berechnen Sie Prognosemodelle. Das sind oft lineare Modelle, also etwa lm(av ~ uv). Es empfiehlt sich, mehrere Modelle zu berechnen und zu schauen, welches dieser Kandidaten am besten ist. Die Güte eines Prognosemodells bemisst sich letztlich nur an der Präzision der Vorhersage neuer Daten, also des Test-Datensatzes. Wie gut Ihre Vorhersagen also wirklich sind, erfahren Sie erst mit der Notenbekanntgabe. Allerdings können Sie die Trainingsdaten nutzen, um die Güte Ihrere Modell abzuschätzen.

2.7.5 Vorhersagen

Schließlich entscheiden Sie sich für einen Modellkandidaten. Diesen Modellkandidaten nehmen Sie her, um die (Ihenn unbekannten) Werte der AV (Zielvariablen) vorherzusagen. Diese Vorhersagen - zusammen mit der ID für jede Vorhersagen - speichern Sie als (reguläre) CSV-Datei ab und reichen Sie als Ihre Prüfungsleistung ein, zusammen mit Ihrer Analysedatei.

2.8.1 Tipps für eine gute Prognose

• Schauen Sie in die Literatur.

• Evtl. kann eine Datenvorverarbeitung (Variablentransformation, z.B. \(\log()\) oder die Elimination von Ausreißern) helfen.

• Überlegen Sie sich Kriterien zur Modell- und/ oder Variablenauswahl. Auch hierfür gibt es Algorithmen und R-Funktionen.

• Vermeiden Sie Über-Anpassung (Overfitting).

• Vermeiden Sie viele fehlende Werte bei Ihrer Prognose. Fehlende Werte werden bei der Benotung mit dem Mittelwert (der vorhandenen Prognosewerte Ihrer Einreichung) aufgefüllt.

• Arbeiten Sie die bereitgestellten Fallstudien durch. Wenn Sie mehr tun möchten, finden Sie im Internet eine Fülle von weiteren Fallstudien.

2.8.2 Tipps zur Datenverarbeitung

• Ein “deutsches” Excel kann Standard-CSV-Dateien nicht ohne Weiteres lesen. Online-Dienste wie Google Sheets können dies allerdings.

2.8.3 Tipps zum Aufbau des Analyseskripts

• Zu Beginn des Skripts sollten alle verwendeten R-Pakete mittels library() gestartet werden.

• Zu Beginn des Skripts sollten die Daten von der vom Dozenten bereitgestellten URL importiert werden (nicht von der eigenen Festplatte, da das Skript sonst bei Dritten, wie Ihrem Prüfer, nicht lauffähig ist).

2.8.4 Sonstiges

• Legen Sie regelmäßig Sicherheitskopien Ihrer Arbeit an (ggf. auf einem anderen Datenträger).

• Achten Sie darauf, dass Sie nicht durcheinander kommen, in welcher Datei der aktuelle Stand Ihrer Arbeit liegt.

2.9.1 Kriterien

• Es gibt drei Bewertungskriterien:

  • Formalia: u.a. Reproduzierbarkeit der Analyse, Lesbarkeit der Syntax, Übersichtlichkeit der Analyse.

  • Methode: u.a. methodischer Anspruch und Korrektheit in der Explorativen Datenanalyse, Datenvorverarbeitung, Variablenauswahl und Modellierungsmethode.

  • Inhalt: Vorhersagegüte.

• Das zentrale Bewertungskriterium ist Inhalt; die übrigen beiden Kriterien fließen nur bei besonders guter oder schlechter Leistung in die Gesamtnote ein.

• Die quantitative Datenanalyse in Durchführung und Interpretation ist der Schwerpunkt dieser Arbeit. Zufälliges identisches Vorgehen, z.B. im R Code, ist sehr unwahrscheinlich und kann als Plagiat bewertet werden.

• Die Gesamtnote muss sich nicht als arithmetischer Mittelwert der Teilnoten ergeben.

• Es werden keine Teilnoten vergeben, sondern nur eine Gesamtnote wird vergeben.

• Es werden keine Hinweise vergeben, stattdessen gibt es einen Überblick an typischen Fehlern.

• Es wird keine Musterlösugn veröffentlicht, um nachfolgende Kohorten nicht zu bevorteilen bzw. die aktuelle Kohorte nicht zu benachteiligen.

2.9.2 Kennzahl der Modellgüte

Die Güte der Vorhersage wird anhand des mittleren Absolutfehlers (mae) bemessen:

\[\text{mae} = \frac{1}{n} \sum_{i=1}^n|(y_i - \hat{y}_i)|\]

2.9.3 Notenstufen

Zur Vorhersagegüte: Die Vorhersagegüte eines einfachen Minimalmodells entspricht einer \(4,0\), die eines Referenzmodells des Dozenten einer \(2,0\).

Ihre Bewertung erfolgt entsprechend Ihrer Vorhersagegüte, d.h., sind Sie besser als das Referenzmodell erhalten Sie hier in diesem Teilaspekt eine bessere Note als \(2,0\)!

2.9.4 Bewertungsprozess

Der Gutachter legt im Nachgang der Prüfung alle Teilnehmis ihre jeweilige Wert der Kennzahl der Modellgüte offen. Außerdem werden die vorherzusagenden Daten veröffentlicht sowie die Grenzwerte für jede Notenstufe. Auf dieser Basis ist es allen Teilnehmis möglich, die Korrektheit Ihrer Note zu überprüfen.

2.12.1 Wo finde ich Beispiele und Vorlagen?

Im Rahmen des Unterrichts wurden mehrere Fallstudien erarbeitet bzw. bereitgestellt, diese dienen Ihnen als ideale Vorlage.

Eine Beispiel-Modellierung finden Sie in der Datei Beispielanalyse-Prognose-Wettbewerb.Rmd. Eine beispielhafte Vorlage (Template), die Sie als Richtschnur nutzen können, ist mit der Datei Template-Vorhersagemodellierung.Rmd hier bereitgestellt.

Im Internet finden sich viele Fallstudien, von denen Sie sich inspirieren lassen können.

2.12.2 Materialsammlung

In [diesem Ordner]((https://github.com/sebastiansauer/Lehre/tree/main/Hinweise/Prognosewettbewerb) finden Sie eine Materialsammlung zum Prognosewettbewerb.

2.12.3 Videos

Diese Playlist beinhaltet Videos, die die Rahmenbedingungen der Prüfungsleistung vorstellt.

3.1.1 Vorbereitung

• Lesen Sie die Hinweise zum Modul.
• Installieren (oder Updaten) Sie die für dieses Modul angegeben Software.
• Lesen Sie die Literatur.

3.1.2 Lernziele

• Sie können erläutern, was man unter statistischem Lernen versteht.
• Sie wissen, war Overfitting ist, wie es entsteht, und wie es vermieden werden kann.
• Sie kennen verschiedenen Arten von statistischem Lernen und können Algorithmen zu diesen Arten zuordnen.

3.1.3 Literatur

• Rhys, Kap. 1
• evtl. Sauer, Kap. 15

3.1.4 Hinweise

• Bitte beachten Sie die Hinweise zum Präsenzunterricht und der Streamingoption.
• Bitte stellen Sie sicher, dass Sie einen einsatzbereiten Computer haben und dass die angegebene Software (in aktueller Version) läuft.

3.3.1 Rule-based

Klassische (ältere) KI implementiert Regeln “hartverdrahtet” in ein Computersystem. Nutzer füttern Daten in dieses System. Das System leitet dann daraus Antworten ab.

Regeln kann man prototypisch mit Wenn-Dann-Abfragen darstellen:

lernzeit <- c(0, 10, 10, 20)
schlauer_nebensitzer <- c(FALSE, FALSE, TRUE, TRUE)

for (i in 1:4) {
  if (lernzeit[i] > 10) {
    print("bestanden!")
  } else {
    if (schlauer_nebensitzer[i] == TRUE) {
      print("bestanden!")
    } else print("Durchgefallen!")
  }
}

## [1] "Durchgefallen!"
## [1] "Durchgefallen!"
## [1] "bestanden!"
## [1] "bestanden!"

Sicherlich könnte man das schlauer programmieren, vielleicht so:

d <- 
  tibble(
  lernzeit = c(0, 10, 10, 20),
  schlauer_nebensitzer = c(FALSE, FALSE, TRUE, TRUE)
)

d %>% 
  mutate(bestanden = ifelse(lernzeit > 10 | schlauer_nebensitzer == TRUE, TRUE, FALSE))

## # A tibble: 4 × 3
##   lernzeit schlauer_nebensitzer bestanden
##      <dbl> <lgl>                <lgl>    
## 1        0 FALSE                FALSE    
## 2       10 FALSE                FALSE    
## 3       10 TRUE                 TRUE     
## 4       20 TRUE                 TRUE

3.3.2 Data-based

ML hat zum Ziel, Regeln aus den Daten zu lernen. Man füttert Daten und Antworten in das System, das System gibt Regeln zurück.

James et al. (2021) definieren ML so: Nehmen wir an, wir haben die abhängige Variable \(Y\) und \(p\) Prädiktoren, \(X_1,X_2, \ldots, X_p\). Weiter nehmen wir an, die Beziehung zwischen \(Y\) und \(X = (X_1, X_2, \ldots, X_p)\) kann durch eine Funktion \(f\) beschrieben werden. Das kann man so darstellen:

\[Y = f(X) + \epsilon\]

ML kann man auffassen als eine Menge an Verfahren, um \(f\) zu schätzen.

Ein Beispiel ist in Abb. 3.3 gezeigt (James et al. 2021).

Abbildung 3.3: Vorhersage des Einkommens durch Ausbildungsjahre

Natürlich kann \(X\) mehr als eine Variable beinhalten, vgl. Abb. 3.4 (James et al. 2021).

Abbildung 3.4: Vorhersage des Einkommens als Funktion von Ausbildungsjahren und Dienstjahren

Anders gesagt: traditionelle KI-Systeme werden mit Daten und Regeln gefüttert und liefern Antworten. ML-Systeme werden mit Daten und Antworten gefüttert und liefern Regeln zurück, vgl. Abb. 3.5.

3.4.1 Modell

Ein Modell, s. Abb. 3.6 (Spurzem 2017)!

Abbildung 3.6: Ein Modell-Auto

Wie man sieht, ist ein Modell eine vereinfachte Repräsentation eines Gegenstands.

Der Gegenstand definiert (gestaltet) das Modell. Das Modell ist eine Vereinfachung des Gegenstands, vgl. Abb. 3.7.

Abbildung 3.7: Gegenstand und Modell

Im maschinellen Lernen meint ein Modell, praktisch gesehen, die Regeln, die aus den Daten gelernt wurden.

3.4.2 Beispiel für einen ML-Algorithmus

Unter einem ML-Algorithmus versteht man das (mathematische oder statistische) Verfahren, anhand dessen die Beziehung zwischen \(X\) und \(Y\) “gelernt” wird. Bei Rhys (2020) (S. 9) findet sich dazu ein Beispiel, das kurz zusammengefasst etwa so lautet:

Beispiel eines Regressionsalgorithmus

1. Setze Gerade in die Daten mit \(b_0 = \hat{y}, b_1 = 0\)
2. Berechne \(MSS = \sum (y_i - \hat{y_i})^2\)
3. “Drehe” die Gerade ein bisschen, d.h. erhöhe \(b_1^{neu} = b_1^{alt} + 0.1\)
4. Wiederhole 2-3 solange, bis \(MSS < \text{Zielwert}\)

Diesen Algorithmus kann man “von Hand” z.B. mit dieser App durchspielen.

3.5.1 Geleitetes Lernen

Die zwei Phasen des geleiteten Lernens sind in Abb. 3.9 dargestellt.

ggplot(mtcars) +
  aes(x = hp, y = mpg) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()

3.5.2 Ungeleitetes Lernen

Die zwei Phasen des ungeleiteten Lernens sind in Abb. 3.10 dargestellt.

Ungeleitetes Lernen kann man wiederum in zwei Arten unterteilen, vgl. Abb. 3.11:

1. Fallreduzierendes Modellieren (Clustering)
2. Dimensionsreduzierendes Modellieren (z.B. Faktorenanalyse)

Abbildung 3.11: Zwei Arten von ungeleitetem Modellieren

4.1.1 Vorbereitung

• Lesen Sie die Literatur.

4.1.2 Lernziele

• Sie können Funktionen, in R schreiben.
• Sie können Datensätze vom Lang- und Breit-Format wechseln.
• Sie können Wiederholungsstrukturen wie Mapping-Funktionen anwenden.
• Sie können eine dplyr-Funktion auf mehrere Spalten gleichzeitig anwenden.

4.1.3 Literatur

• Rhys, Kap. 2
• MODAR, Kap. 5

4.2.1 Überblick

In R ist praktisch alles ein Objekt. Ein Objekt meint ein im Computerspeicher repräsentiertes Ding, etwa eine Tabelle.

Vektoren und Dataframes (Tibbles) sind die vielleicht gängigsten Objektarten in R (vgl. Abb. 4.1, aus Sauer (2019)).

Abbildung 4.1: Zentrale Objektarten in R

Es gibt in R keine (Objekte für) Skalare (einzelne Zahlen). Stattdessen nutzt R Vektoren der Länge 1.

Ein nützliches Schema stammt aus Wickham and Grolemund (2016), s. Abb. 4.2.

Abbildung 4.2: Objektarten hierarchisch gegliedert

4.2.2 Taxonomie

Unter homogenen Objektiven verstehen wir Datenstrukturen, die nur eine Art von Daten (wie Text oder Ganze Zahlen) fassen. Sonstige Objekte nennen wir heterogen.

• Homogene Objekte
  • Vektoren
  • Matrizen
• Heterogen
  • Liste
  • Dataframes (Tibbles)

4.2.2.1 Vektoren

Vektoren sind insofern zentral in R, als dass die übrigen Datenstrukturen auf ihnen aufbauen, vgl. Abb. 4.3 aus Sauer (2019).

Reine (atomare) Vektoren in R sind eine geordnete Liste von Daten eines Typs.

Abbildung 4.3: Vektoren stehen im Zentrum der Datenstrukturen in R

ein_vektor <- c(1, 2, 3)
noch_ein_vektor <- c("A", "B", "C")
logischer_vektor <- c(TRUE, FALSE, TRUE)

Mit str() kann man sich die Struktur eines Objektsausgeben lassen:

str(ein_vektor)

##  num [1:3] 1 2 3

str(noch_ein_vektor)

##  chr [1:3] "A" "B" "C"

str(logischer_vektor)

##  logi [1:3] TRUE FALSE TRUE

Vektoren können von folgenden Typen sein:

• Kommazahlen ( double) genannt
• Ganzzahlig (integer, auch mit L für Long abgekürzt)
• Text (´character`, String)
• logische Ausdrücke (logical oder lgl) mit TRUE oder FALSE

Kommazahlen und Ganze Zahlen zusammen bilden den Typ numeric (numerisch) in R.

knitr::opts_chunk$set(echo = TRUE)

Den Typ eines Vektors kann man mit typeof() ausgeben lassen:

typeof(ein_vektor)

## [1] "double"

sex <- factor(c("Mann", "Frau", "Frau"))

str(sex)

##  Factor w/ 2 levels "Frau","Mann": 2 1 1

eine_liste <- list(titel = "Einführung",
                   woche = 1,
                   datum = c("2022-03-14", "2202-03-21"),
                   lernziele = c("dies", "jenes", "und noch mehr"),
                   lehre = c(TRUE, TRUE, TRUE)
                   )
str(eine_liste)

## List of 5
##  $ titel    : chr "Einführung"
##  $ woche    : num 1
##  $ datum    : chr [1:2] "2022-03-14" "2202-03-21"
##  $ lernziele: chr [1:3] "dies" "jenes" "und noch mehr"
##  $ lehre    : logi [1:3] TRUE TRUE TRUE

library(tidyverse)

studentis <-
  tibble(
    name = c("Anna", "Berta"),
    motivation = c(10, 20),
    noten = c(1.3, 1.7)
  )
str(studentis)

## tibble [2 × 3] (S3: tbl_df/tbl/data.frame)
##  $ name      : chr [1:2] "Anna" "Berta"
##  $ motivation: num [1:2] 10 20
##  $ noten     : num [1:2] 1.3 1.7

4.2.3 Indizieren

Einen Teil eines Objekts auszulesen, bezeichnen wir als Indizieren.

str(ein_vektor)

##  num [1:3] 1 2 3

ein_vektor[1]

## [1] 1

ein_vektor[c(1,2)]

## [1] 1 2

ein_vektor[1,2]

## Error in ein_vektor[1, 2]: incorrect number of dimensions

ein_vektor[[2]]

## [1] 2

is.vector(eine_liste)

## [1] TRUE

ein_vektor[4]

## [1] NA

eine_liste %>% str()

## List of 5
##  $ titel    : chr "Einführung"
##  $ woche    : num 1
##  $ datum    : chr [1:2] "2022-03-14" "2202-03-21"
##  $ lernziele: chr [1:3] "dies" "jenes" "und noch mehr"
##  $ lehre    : logi [1:3] TRUE TRUE TRUE

eine_liste[1]

## $titel
## [1] "Einführung"

eine_liste[2]

## $woche
## [1] 1

eine_liste[[1]]

## [1] "Einführung"

eine_liste[["titel"]]

## [1] "Einführung"

eine_liste$titel

## [1] "Einführung"

eine_liste[[4]][[1]] 

## [1] "dies"

pluck(eine_liste, 4)

## [1] "dies"          "jenes"         "und noch mehr"

pluck(eine_liste, 4, 1)

## [1] "dies"

eine_liste %>% length()

## [1] 5

eine_liste[[6]]

## Error in eine_liste[[6]]: subscript out of bounds

eine_liste[6]

## $<NA>
## NULL

studentis[1]

## # A tibble: 2 × 1
##   name 
##   <chr>
## 1 Anna 
## 2 Berta

studentis["name"]

## # A tibble: 2 × 1
##   name 
##   <chr>
## 1 Anna 
## 2 Berta

studentis[["name"]]

## [1] "Anna"  "Berta"

studentis$name

## [1] "Anna"  "Berta"

4.2.4 Weiterführende Hinweise

• Tutorial zum Themen Indizieren von Listen von Jenny BC.

4.2.5 Indizieren mit dem Tidyverse

Natürlich kann man auch die Tidyverse-Verben zum Indizieren verwenden. Das bietet sich an, wenn zwei Bedingungen erfüllt sind:

1. Wenn man einen Tibble als Input und als Output hat
2. Wenn man nicht programmieren möchte

4.5.1 across()

Handelt es sich um Spalten von Tibbles, dann bietet sich die Funktion across(.col, .fns) an. across wendet eine oder mehrere Funktionen (mit .fns bezeichnet) auf die Spalten .col an.

Das erklärt sich am besten mit einem Beispiel:

Natürlich hätte man in diesem Fall auch anders vorgehen können:

mtcars %>% 
  summarise(across(.cols = everything(),
                   .fns = mean))

##        mpg    cyl     disp       hp     drat      wt     qsec     vs      am
## 1 20.09062 6.1875 230.7219 146.6875 3.596563 3.21725 17.84875 0.4375 0.40625
##     gear   carb
## 1 3.6875 2.8125

Möchte man der Funktion .fns Parameter übergeben, so nutzt man diese Syntax (“Purrr-Lambda”):

mtcars %>% 
  summarise(across(.cols = everything(),
                   .fns = ~ mean(., na.rm = TRUE)))

##        mpg    cyl     disp       hp     drat      wt     qsec     vs      am
## 1 20.09062 6.1875 230.7219 146.6875 3.596563 3.21725 17.84875 0.4375 0.40625
##     gear   carb
## 1 3.6875 2.8125

Hier findet sich ein guter Überblick zu across().

4.5.2 map()

map() ist eine Funktion aus dem R-Paket purrr und Teil des Tidyverse.

map(x, f) wenden die Funktion f auf jedes Element von x an. Ist x ein Tibble, so wird f demnach auf jede Spalte von x angewendet (“zugeordnet”, daher map), vgl. Abb. 4.6 aus Sauer (2019).

Abbildung 4.6: Sinnbild für map

Hier ein Beispiel-Code:

data(mtcars)

mtcars <- mtcars %>% select(1:3)  # nur die ersten 3 Spalten

map(mtcars, mean)

## $mpg
## [1] 20.09062
## 
## $cyl
## [1] 6.1875
## 
## $disp
## [1] 230.7219

Möchte man der gemappten Funktion Parameter übergeben, nutzt man wieder die “Kringel-Schreibweise”:

map(mtcars, ~ mean(., na.rm = TRUE))

## $mpg
## [1] 20.09062
## 
## $cyl
## [1] 6.1875
## 
## $disp
## [1] 230.7219

4.5.3 Weiterführende Hinweise

Weiteres zu map() findet sich z.B. in Wickham and Grolemund (2016), Kapitel 21.5 oder in Sauer (2019), Kap. 28.2.

Tutorial zu map() von Jenny BC.

4.6.1 Wozu Listenspalten?

Listenspalten sind immer dann sinnvoll, wenn eine einfache Tabelle nicht komplex genug für unsere Daten ist.

Zwei Fälle stechen dabei ins Auge:

1. Unsere Datenstruktur ist nicht rechteckig
2. In einer Zelle der Tabelle soll mehr als ein einzelner Wert stehen: vielleicht ein Vektor, eine Liste oder eine Tabelle

Der erstere Fall (nicht reckeckig) ließe sich noch einfach lösen, in dem man mit NA auffüllt.

Der zweite Fall verlangt schlichtweg nach komplexeren Datenstrukturen.

Kap. 25.3 aus Wickham and Grolemund (2016) bietet einen guten Einstieg in das Konzept von Listenspalten (list-columns) in R.

4.6.2 Beispiele für Listenspalten

library(tidymodels)
mtcars_cv <-
  vfold_cv(mtcars, v = 3)

mtcars_cv

## #  3-fold cross-validation 
## # A tibble: 3 × 2
##   splits          id   
##   <list>          <chr>
## 1 <split [21/11]> Fold1
## 2 <split [21/11]> Fold2
## 3 <split [22/10]> Fold3

pos11 <- mtcars_cv[[1]][[1]]
pos11

## <Analysis/Assess/Total>
## <21/11/32>

str(pos11)

## List of 4
##  $ data  :'data.frame':  32 obs. of  3 variables:
##   ..$ mpg : num [1:32] 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##   ..$ cyl : num [1:32] 6 6 4 6 8 6 8 4 4 6 ...
##   ..$ disp: num [1:32] 160 160 108 258 360 ...
##  $ in_id : int [1:21] 1 2 6 8 9 11 12 15 17 18 ...
##  $ out_id: logi NA
##  $ id    : tibble [1 × 1] (S3: tbl_df/tbl/data.frame)
##   ..$ id: chr "Fold1"
##  - attr(*, "class")= chr [1:2] "vfold_split" "rsplit"

pos11[["data"]][[1]]

##  [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 10.4
## [16] 10.4 14.7 32.4 30.4 33.9 21.5 15.5 15.2 13.3 19.2 27.3 26.0 30.4 15.8 19.7
## [31] 15.0 21.4

pos11[["data"]][1] %>% 
  head()

##                    mpg
## Mazda RX4         21.0
## Mazda RX4 Wag     21.0
## Datsun 710        22.8
## Hornet 4 Drive    21.4
## Hornet Sportabout 18.7
## Valiant           18.1

pluck(pos11, "data", 1, 1)

## [1] 21

tibble::tribble(
  ~Nr,   ~Titel,   ~Literatur,   ~Aufgaben,
   1L, "Titel1", "Literatur1", "Aufgaben1",
   2L, "Titel2", "Literatur2", "Aufgaben2",
   3L, "Titel3", "Literatur3", "Aufgaben3"
  ) %>% 
  gt::gt()

source("https://raw.githubusercontent.com/sebastiansauer/Lehre/main/R-Code/render-course-sections.R")

course_dates_file <- "https://raw.githubusercontent.com/sebastiansauer/datascience1/main/course-dates.yaml"

content_file <- "https://raw.githubusercontent.com/sebastiansauer/datascience1/main/_modul-ueberblick.yaml"

mastertable <- build_master_course_table(
  course_dates_file =  course_dates_file,
  content_file = content_file)

mastertable[["Vorbereitung"]][[1]]

## [1] "Lesen Sie die Hinweise zum Modul."                                       
## [2] "Installieren (oder Updaten) Sie die für dieses Modul angegeben Software."
## [3] "Lesen Sie die Literatur."

mastertable[["Vorbereitung"]][[1]][2]

## [1] "Installieren (oder Updaten) Sie die für dieses Modul angegeben Software."

mastertable["Vorbereitung"] %>% class()

## [1] "tbl_df"     "tbl"        "data.frame"

mastertable[["Vorbereitung"]] %>% class()

## [1] "list"

mastertable[["Vorbereitung"]][[1]]

## [1] "Lesen Sie die Hinweise zum Modul."                                       
## [2] "Installieren (oder Updaten) Sie die für dieses Modul angegeben Software."
## [3] "Lesen Sie die Literatur."

mastertable[["Vorbereitung"]][1]

## [[1]]
## [1] "Lesen Sie die Hinweise zum Modul."                                       
## [2] "Installieren (oder Updaten) Sie die für dieses Modul angegeben Software."
## [3] "Lesen Sie die Literatur."

mastertable[["Vorbereitung"]][1] %>% pluck(2)

## NULL

mastertable[["Vorbereitung"]][1][2]

## [[1]]
## NULL

mastertable[["Vorbereitung"]][1][[2]]

## Error in mastertable[["Vorbereitung"]][1][[2]]: subscript out of bounds

mastertable[["Vorbereitung"]][[1]][2]

## [1] "Installieren (oder Updaten) Sie die für dieses Modul angegeben Software."

4.6.3 Programmieren mit dem Tidyverse

Das Programmieren mit dem Tidyvers ist nicht ganz einfach und hier nicht näher ausgeführt. Eine Einführung findet sich z.B.

• Tidyeval in fünf Minuten (Video)
• In Kapiteln 17-21 in Advanced R, 2nd Ed
• Ein Überblicksdiagramm findet sich hier Quelle.

5.1.1 Vorbereitung

• Lesen Sie TMWR, Kapitel 1
• Lesen Sie übrige Literatur zu diesem Thema

5.1.2 Lernziele

• Sie sind in der Lage, Regressionsmodelle mit dem tidymodels-Ansatz zu spezifizieren

5.1.3 Literatur

• TMWR, Kap. 1, 5, 6, 7, 8, 9

5.4.1 Modelle spezifizieren

Ein (statistisches) Modell wird in Tidymodels mit drei Elementen spezifiziert, vgl. Abb. 5.1.

Die Definition eines Modells in tidymodels folgt diesen Ideen:

1. Das Modell sollte unabhängig von den Daten spezifiziert sein
2. Das Modell sollte unabhängig von den Variablen (AV, UVs) spezifiziert sein
3. Das Modell sollte unabhängig von etwaiger Vorverarbeitung (z.B. z-Transformation) spezifiziert sein

Da bei einer linearen Regression nur der Modus “Regression” möglich ist, muss der Modus in diesem Fall nicht angegeben werden. Tidymodels erkennt das automatisch.

lm_model <-   
  linear_reg() %>%   # Algorithmus, Modelltyp
  set_engine("lm")  # Implementierung
  # Modus hier nicht nötig, da lineare Modelle immer numerisch klassifizieren

5.4.2 Modelle berechnen

Nach Rhys (2020) ist ein Modell sogar erst ein Modell, wenn die Koeffizienten berechnet sind. Tidymodels kennt diese Unterscheidung nicht. Stattdessen spricht man in Tidymodels von einem “gefitteten” Modell, sobald es berechnet ist. Ähnlich fancy könnte man von einem “instantiierten” Modell sprechen.

Für das Beispiel der einfachen linearen Regression heißt das, das Modell ist gefittet, sobald die Steigung und der Achsenabschnitt (sowie die Residualstreuung) berechnet sind.

lm_form_fit <- 
  lm_model %>% 
  fit(Sale_Price ~ Longitude + Latitude, data = ames_train)

5.4.3 Vorhersagen

Im maschinellen Lernen ist man primär an den Vorhersagen interessiert, häufig nur an Punktschätzungen. Schauen wir uns also zunächst diese an.

Vorhersagen bekommt man recht einfach mit der predict() Methode:

predict(lm_form_fit, new_data = ames_test) %>% 
  head()

## # A tibble: 6 × 1
##   .pred
##   <dbl>
## 1  5.23
## 2  5.29
## 3  5.29
## 4  5.28
## 5  5.24
## 6  5.24

Die Syntax lautet predict(modell, daten_zum_vorhersagen).

5.4.4 Vorhersagen im Train-Datensatz

Vorhersagen im Train-Datensatz machen keinen Sinn, da sie nicht gegen Overfitting geschützt sind und daher deutlich zu optimistisch sein können.

Bei einer linearen Regression ist diese Gefahr nicht so hoch, aber bei anderen, flexibleren Modellen, ist diese Gefahr absurd groß.

5.4.5 Modellkoeffizienten im Train-Datensatz

Gibt man den Namen des Modellobjekts ein, so wird ein Überblick an relevanten Modellergebnissen am Bildschirm gedruckt:

lm_form_fit

## parsnip model object
## 
## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -312.289       -2.094        2.888

Innerhalb des Ergebnisobjekts findet sich eine Liste namens fit, in der die Koeffizienten (der “Fit”) abgelegt sind:

lm_form_fit %>% pluck("fit")

## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -312.289       -2.094        2.888

Zum Herausholen dieser Infos kann man auch die Funktion extract_fit_engine() verwenden:

lm_fit <-
  lm_form_fit %>% 
  extract_fit_engine()

lm_fit

## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -312.289       -2.094        2.888

Das extrahierte Objekt ist, in diesem Fall, das typische lm() Objekt. Entsprechend kann man daruaf coef() oder summary() anwenden.

coef(lm_fit)

## (Intercept)   Longitude    Latitude 
## -312.288677   -2.094429    2.887655

summary(lm_fit)

## 
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02650 -0.09839 -0.01662  0.09653  0.57653 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -312.2887    14.7509  -21.17   <2e-16 ***
## Longitude     -2.0944     0.1324  -15.82   <2e-16 ***
## Latitude       2.8877     0.1815   15.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1623 on 2339 degrees of freedom
## Multiple R-squared:  0.1728, Adjusted R-squared:  0.1721 
## F-statistic: 244.4 on 2 and 2339 DF,  p-value: < 2.2e-16

Schicker sind die Pendant-Befehle aus broom, die jeweils einen Tibble zuückliefern:

library(broom)
tidy(lm_fit) # Koeffizienten

## # A tibble: 3 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  -312.      14.8       -21.2 3.71e-91
## 2 Longitude      -2.09     0.132     -15.8 1.25e-53
## 3 Latitude        2.89     0.181      15.9 3.36e-54

glance(lm_fit) # Modellgüte

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic  p.value    df logLik    AIC    BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl>  <dbl>  <dbl>
## 1     0.173         0.172 0.162      244. 4.20e-97     2   937. -1866. -1843.
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

5.4.6 Parsnip RStudio add-in

Mit dem Add-in von Parsnip kann man sich eine Modellspezifikation per Klick ausgeben lassen. Nett!

parsnip_addin()

5.5.1 Konzept des Workflows in Tidymodels

5.5.2 Einfaches Beispiel

Wir initialisieren einen Workflow, verzichten auf Vorverarbeitung und fügen ein Modell hinzu:

lm_workflow <- 
  workflow() %>%  # init
  add_model(lm_model) %>%   # Modell hinzufügen
  add_formula(Sale_Price ~ Longitude + Latitude)  # Modellformel hinzufügen

Werfen wir einen Blick in das Workflow-Objekt:

lm_workflow

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

Wie man sieht, gehört die Modellformel (y ~ x) zur Vorverarbeitung aus Sicht von Tidymodels.

Was war nochmal im Objekt lm_model enthalten?

lm_model

## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

Jetzt können wir das Modell berechnen (fitten):

lm_fit <- 
  lm_workflow %>%
  fit(ames_train)

Natürlich kann man synonym auch schreiben:

lm_fit <- fit(lm_wflow, ames_train)

Schauen wir uns das Ergebnis an:

lm_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:
## stats::lm(formula = ..y ~ ., data = data)
## 
## Coefficients:
## (Intercept)    Longitude     Latitude  
##    -312.289       -2.094        2.888

5.5.3 Vorhersage mit einem Workflow

Die Vorhersage mit einem Tidymodels-Workflow ist einerseits komfortabel, da man einfach sagen kann:

“Nimm die richtigen Koeffizienten des Modells aus dem Train-Set und wende sie auf das Test-Sample an. Berechne mir die Vorhersagen und die Modellgüte.”

So sieht das aus:

final_lm_res <- last_fit(lm_workflow, ames_split)
final_lm_res

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits             id               .metrics .notes   .predictions .workflow 
##   <list>             <chr>            <list>   <list>   <list>       <list>    
## 1 <split [2342/588]> train/test split <tibble> <tibble> <tibble>     <workflow>

Anderseits wird auch ein recht komplexes Objekt zurückgeliefert, das man erst mal durchschauen muss.

Wie man sieht, gibt es mehrere Listenspalten. Besonders interessant erscheinen natürlich die Listenspalten .metrics und .predictions.

Schauen wir uns die Vorhersagen an.

lm_preds <- final_lm_res %>% pluck(".predictions", 1)

Es gibt auch eine Funktion, die obige Zeile vereinfacht (also synonym ist):

lm_preds <- collect_predictions(final_lm_res)
lm_preds %>% slice_head(n = 5)

## # A tibble: 5 × 5
##   id               .pred  .row Sale_Price .config             
##   <chr>            <dbl> <int>      <dbl> <chr>               
## 1 train/test split  5.23     2       5.02 Preprocessor1_Model1
## 2 train/test split  5.29     5       5.28 Preprocessor1_Model1
## 3 train/test split  5.29     6       5.29 Preprocessor1_Model1
## 4 train/test split  5.28     9       5.37 Preprocessor1_Model1
## 5 train/test split  5.24    27       5.10 Preprocessor1_Model1

5.5.4 Modellgüte

Dieser Abschnitt bezieht sich auf Kapitel 9 in Silge and Kuhn (2022).

Die Vorhersagen bilden die Basis für die Modellgüte (“Metriken”), die schon fertig berechnet im Objekt final_lm_res liegen und mit collect_metrics herausgenommen werden können:

lm_metrics <- collect_metrics(final_lm_res)

Man kann auch angeben, welche Metriken der Modellgüte man bekommen möchte:

ames_metrics <- metric_set(rmse, rsq, mae)

ames_metrics(data = lm_preds, 
             truth = Sale_Price, 
             estimate = .pred)

## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       0.156
## 2 rsq     standard       0.171
## 3 mae     standard       0.121

5.5.5 Vorhersage von Hand

Man kann sich die Metriken auch von Hand ausgeben lassen, wenn man direktere Kontrolle haben möchte als mit last_fit und collect_metrics.

ames_test_small <- ames_test %>% slice(1:5)
predict(lm_form_fit, new_data = ames_test_small)

## # A tibble: 5 × 1
##   .pred
##   <dbl>
## 1  5.23
## 2  5.29
## 3  5.29
## 4  5.28
## 5  5.24

Jetzt binden wir die Spalten zusammen, also die “Wahrheit” (\(y\)) und die Vorhersagen:

ames_test_small2 <- 
  ames_test_small %>% 
  select(Sale_Price) %>% 
  bind_cols(predict(lm_form_fit, ames_test_small)) %>% 
  # Add 95% prediction intervals to the results:
  bind_cols(predict(lm_form_fit, ames_test_small, type = "pred_int")) 

rsq(ames_test_small2, 
   truth = Sale_Price,
   estimate = .pred
   )

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rsq     standard       0.832

Andere Koeffizienten der Modellgüte können mit rmse oder mae abgerufen werden.

5.6.1 Was ist Rezept und wozu ist es gut?

So könnte ein typischer Aufruf von lm() aussehen:

lm(Sale_Price ~ Neighborhood + log10(Gr_Liv_Area) + Year_Built + Bldg_Type, 
   data = ames)

Neben dem Fitten des Modells besorgt die Formel-Schreibweise noch einige zusätzliche nützliche Vorarbeitung:

1. Definition von AV und AV
2. Log-Transformation von Gr_Liv_Area
3. Transformation der nominalen Variablen in Dummy-Variablen

Das ist schön und nütlich, hat aber auch Nachteile:

1. Das Modell wird nicht nur spezifiziert, sondern auch gleich berechnet. Das ist unpraktisch, weil man die Modellformel vielleicht in anderen Modell wiederverwenden möchte. Außerdem kann das Berechnen lange dauern.
2. Die Schritte sind ineinander vermengt, so dass man nicht einfach und übersichtlich die einzelnen Schritte bearbeiten kann.

Praktischer wäre also, die Schritte der Vorverarbeitung zu ent-flechten. Das geht mit einem “Rezept” aus Tidmoodels:

simple_ames <- 
  recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area + Year_Built + Bldg_Type,
         data = ames_train) %>%
  step_log(Gr_Liv_Area, base = 10) %>% 
  step_dummy(all_nominal_predictors())
simple_ames

## Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor          4
## 
## Operations:
## 
## Log transformation on Gr_Liv_Area
## Dummy variables from all_nominal_predictors()

5.6.2 Workflows mit Rezepten

Jetzt definieren wir den Workflow nicht nur mit einer Modellformel, sondern mit einem Rezept:

lm_workflow <-
  workflow() %>% 
  add_model(lm_model) %>% 
  add_recipe(simple_ames)

Sonst hat sich nichts geändert.

Wie vorher, können wir jetzt das Modell berechnen.

lm_fit <- fit(lm_workflow, ames_train)

final_lm_res <- last_fit(lm_workflow, ames_split)
final_lm_res

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits             id               .metrics .notes   .predictions .workflow 
##   <list>             <chr>            <list>   <list>   <list>       <list>    
## 1 <split [2342/588]> train/test split <tibble> <tibble> <tibble>     <workflow>
## 
## There were issues with some computations:
## 
##   - Warning(s) x1: prediction from a rank-deficient fit may be misleading
## 
## Use `collect_notes(object)` for more information.

lm_metrics <- collect_metrics(final_lm_res)
lm_metrics

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard      0.0761 Preprocessor1_Model1
## 2 rsq     standard      0.804  Preprocessor1_Model1

5.6.3 Spaltenrollen

Eine praktische Funktion ist es, bestimmte Spalten nicht als Prädiktor, sondern als ID-Variable zu nutzen. Das kann man in Tidymodels komfortabel wie folgt angeben:

ames_recipe <-
  simple_ames %>% 
  update_role(Neighborhood, new_role = "id")

ames_recipe

## Recipe
## 
## Inputs:
## 
##       role #variables
##         id          1
##    outcome          1
##  predictor          3
## 
## Operations:
## 
## Log transformation on Gr_Liv_Area
## Dummy variables from all_nominal_predictors()

5.6.4 Fazit

Mehr zu Rezepten findet sich hier. Ein Überblick zu allen Schritten der Vorverarbeitung findet sich hier.

6.1.1 Lernziele

• Sie sind in der Lage, einfache Klassifikationsmodelle zu spezifizieren mit tidymodels
• Sie können den knn-Algorithmus erläutern
• Sie können den knn-Algorithmus in tidymodels anwenden
• Sie können die Gütemetriken von Klassifikationsmodellen einschätzen

6.1.2 Literatur

• Rhys, Kap. 3
• Timbers et al., Kap. 5

6.6.1 Analog zu Timbers et al.

Eine Anwendung von kNN mit Tidymodels ist in Timbers, Campbell, and Lee (2022), Kap. 5.6, hier beschrieben.

Die Daten aus Timbers, Campbell, and Lee (2022) finden sich in diesem Github-Repo-

Die (z-transformierten) Daten zur Tumorklassifikation können hier bezogen werden.

data_url <- "https://raw.githubusercontent.com/UBC-DSCI/introduction-to-datascience/master/data/wdbc.csv"
cancer <- read_csv(data_url)

Timbers, Campbell, and Lee (2022) verwenden in Kap. 5 auch noch nicht standardisierte Daten, unscales_wdbc.csv, die hier als CSV-Datei heruntergeladen werden können.

cancer_unscales_path <- "https://raw.githubusercontent.com/UBC-DSCI/introduction-to-datascience/master/data/unscaled_wdbc.csv"

unscaled_cancer <- read_csv(cancer_unscales_path) |>
  mutate(Class = as_factor(Class)) |>
  select(Class, Area, Smoothness)
unscaled_cancer

## # A tibble: 569 × 3
##    Class  Area Smoothness
##    <fct> <dbl>      <dbl>
##  1 M     1001      0.118 
##  2 M     1326      0.0847
##  3 M     1203      0.110 
##  4 M      386.     0.142 
##  5 M     1297      0.100 
##  6 M      477.     0.128 
##  7 M     1040      0.0946
##  8 M      578.     0.119 
##  9 M      520.     0.127 
## 10 M      476.     0.119 
## # … with 559 more rows

6.6.2 Rezept definieren

uc_recipe <- recipe(Class ~ ., data = unscaled_cancer)
print(uc_recipe)

## Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor          2

Und jetzt die z-Transformation:

uc_recipe <- 
  uc_recipe |>
  step_scale(all_predictors()) |>
  step_center(all_predictors())

Die Schritte prep() und bake() sparen wir uns, da fit() und predict() das für uns besorgen.

6.6.3 Modell definieren

knn_spec <- nearest_neighbor(weight_func = "rectangular", neighbors = 5) |>
  set_engine("kknn") |>
  set_mode("classification")
knn_spec

## K-Nearest Neighbor Model Specification (classification)
## 
## Main Arguments:
##   neighbors = 5
##   weight_func = rectangular
## 
## Computational engine: kknn

6.6.4 Workflow definieren

knn_fit <- workflow() |>
  add_recipe(uc_recipe) |>
  add_model(knn_spec) |>
  fit(data = unscaled_cancer)

knn_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: nearest_neighbor()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 2 Recipe Steps
## 
## • step_scale()
## • step_center()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:
## kknn::train.kknn(formula = ..y ~ ., data = data, ks = min_rows(5,     data, 5), kernel = ~"rectangular")
## 
## Type of response variable: nominal
## Minimal misclassification: 0.1107206
## Best kernel: rectangular
## Best k: 5

6.6.5 Vorhersagen

new_observation <- tibble(Area = c(500, 1500), Smoothness = c(0.075, 0.1))
prediction <- predict(knn_fit, new_observation)

prediction

## # A tibble: 2 × 1
##   .pred_class
##   <fct>      
## 1 B          
## 2 M

6.7.1 Rezept definieren

cancer_recipe <- recipe(Class ~ Smoothness + Concavity, data = cancer_train) |>
  step_scale(all_predictors()) |>
  step_center(all_predictors())

6.7.2 Modell definieren

knn_spec <- nearest_neighbor(weight_func = "rectangular", neighbors = 3) |>
  set_engine("kknn") |>
  set_mode("classification")

6.7.3 Workflow definieren

knn_fit <- workflow() |>
  add_recipe(cancer_recipe) |>
  add_model(knn_spec) |>
  fit(data = cancer_train)

knn_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: nearest_neighbor()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 2 Recipe Steps
## 
## • step_scale()
## • step_center()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:
## kknn::train.kknn(formula = ..y ~ ., data = data, ks = min_rows(3,     data, 5), kernel = ~"rectangular")
## 
## Type of response variable: nominal
## Minimal misclassification: 0.1197183
## Best kernel: rectangular
## Best k: 3

6.7.4 Vorhersagen

Im Gegensatz zu Timbers, Campbell, and Lee (2022) verwenden wir hier last_fit() und collect_metrics(), da wir dies bereits eingeführt haben und künftig darauf aufbauen werden.

cancer_test_fit <- last_fit(knn_fit, cancer_split)

cancer_test_fit

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits            id               .metrics .notes   .predictions .workflow 
##   <list>            <chr>            <list>   <list>   <list>       <list>    
## 1 <split [426/143]> train/test split <tibble> <tibble> <tibble>     <workflow>

6.7.5 Modellgüte

cancer_test_fit %>% collect_metrics()

## # A tibble: 2 × 4
##   .metric  .estimator .estimate .config             
##   <chr>    <chr>          <dbl> <chr>               
## 1 accuracy binary         0.902 Preprocessor1_Model1
## 2 roc_auc  binary         0.937 Preprocessor1_Model1

Die eigentlichen Predictions stecken in der Listenspalte .predictions im Fit-Objekt:

names(cancer_test_fit)

## [1] "splits"       "id"           ".metrics"     ".notes"       ".predictions"
## [6] ".workflow"

Genau genommen ist .predictions eine Spalte, in der in jeder Zeile (und damit Zelle) eine Tabelle (Tibble) steht. Wir haben nur eine Zeile und wollen das erste Element dieser Spalte herausziehen. Da hilft pluck():

cancer_test_predictions <- 
cancer_test_fit %>% 
  pluck(".predictions", 1)

confusion <- cancer_test_predictions |>
             conf_mat(truth = Class, estimate = .pred_class)

confusion

##           Truth
## Prediction  B  M
##          B 82  6
##          M  8 47

6.7.6 Visualisierung

autoplot(confusion, type = "mosaic")
autoplot(confusion, type = "heatmap") +
  labs(x = "Beobachtung",
       y = "Vorhersage",
       title = "Konfusionsmatrix")

7.1.1 Vorbereitung

• Lesen Sie die Literatur.

7.1.2 Lernziele

• Sie verstehen den Nutzen von Resampling und Tuning im maschinellen Nutzen.
• Sie können Methoden des Resampling und Tunings mit Hilfe von Tidymodels anwenden.

7.1.3 Literatur

• Rhys, Kap. 3
• TMWR, Kap. 10, 12

7.3.1 Datensatz aufteilen

data(ames)

set.seed(4595)
data_split <- initial_split(ames, strata = "Sale_Price")

ames_train <- training(data_split)
ames_test <- testing(data_split)

7.3.2 Rezept, Modell und Workflow definieren

In gewohnter Weise definieren wir den Workflow mit einem kNN-Modell.

ames_rec <-
  recipe(Sale_Price ~ ., data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_other(Neighborhood, threshold = .1)  %>%
  step_dummy(all_nominal()) %>%
  step_zv(all_predictors()) 

knn_model <-
  nearest_neighbor(
    mode = "regression",
  ) %>%
  set_engine("kknn")

ames_wflow <-
  workflow() %>%
  add_recipe(ames_rec) %>%
  add_model(knn_model)

Das kNN-Modell ist noch nicht berechnet, es ist nur ein “Rezept” erstellt:

knn_model

## K-Nearest Neighbor Model Specification (regression)
## 
## Computational engine: kknn

ames_wflow

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: nearest_neighbor()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 4 Recipe Steps
## 
## • step_log()
## • step_other()
## • step_dummy()
## • step_zv()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## K-Nearest Neighbor Model Specification (regression)
## 
## Computational engine: kknn

7.5.1 Einfache v-fache Kreuzvalidierung

Abb. 7.2 illustriert die zufällige Aufteilung von \(n=10\) Fällen der Originalstrichprobe auf eine Train- bzw. Test-Stichpobe. Man spricht von Kreuzvalidierung (cross validation, CV).

In diesem Fall wurden 70% der (\(n=10\)) Fälle der Train-Stichprobe zugewiesen (der Rest der Test-Stichprobe); ein willkürlicher, aber nicht unüblicher Anteil. Diese Aufteilung wurde \(v=3\) Mal vorgenommen, es resultieren drei “Resampling-Stichproben”, die manchmal auch als “Faltungen” bezeichnet werden.

Abbildung 7.2: Resampling: Eine Stichprobe wird mehrfach (hier 3 Mal) zu 70% in ein Train- und zu 30% in die Test-Stichprobe aufgeteilt

Sauer (2019) stellt das Resampling so dar (S. 259), s. Abb. 7.3.

Abbildung 7.3: Kreuzvalidierung, Aufteilung in Train- vs. Testsample

Der Gesamtfehler der Vorhersage wird als Mittelwerte der Vorhersagefehler in den einzelnen Faltungen berechnet.

Warum ist die Vorhersage besser, wenn man mehrere Faltungen, mehrere Schätzungen für \(y\) also, vornimmt?

Der Grund ist das Gesetz der großen Zahl, nachdem sich eine Schätzung in Mittelwert und Variabilität stabilisiert mit steigendem Stichprobenumfang, dem wahren Mittelwert also präziser schätzt. Bei Normalverteilungen klappt das gut, bei randlastigen Verteilungen leider nicht mehr (Taleb 2019).

Häufig werden \(v=10\) Faltungen verwendet, was sich empirisch als guter Kompromiss von Rechenaufwand und Fehlerreduktion herausgestellt hat.

7.5.2 Wiederholte Kreuzvalidierung

Die \(r\)-fach wiederholte Kreuzvalidierung wiederholte die einfache Kreuzvalidierung mehrfach (nämlich \(r=4\) mal), Sauer (2019) stellt das Resampling so dar (S. 259), s. Abb. 7.4.

Abbildung 7.4: Wiederholte Kreuzvalidierung

Die wiederholte Kreuzvalidierung reduziert den Standardfehler der Vorhersagen.

Silge and Kuhn (2022) zeigen die Verringerung des Schätzfehlers als Funktion der \(r\) Wiederholungen dar, s. Abb. 7.5.

Abbildung 7.5: Reduktion des Schätzfehlers als Funktion der r Wiederhoulugen der Kreuzvalidierung

Warum ist die Wiederholung der Kreuzvalidierung nützlich?

Die Kreuvalidierung liefert einen Schätzwert der Modellparameter, die wahren Modellparameter werden also anhand einer Stichprobe von \(n=1\) geschätzt. Mit höherem Stichprobenumfang kann diese Schätzung natürlich präzisiert werden.

Da jede Stichprobenverteilung bei \(n \rightarrow \infty\) normalverteilt ist - ein zentrales Theorem der Statistik, der Zentrale Grenzwertsatz (Central Limit Theorem) - kann man hoffen, dass sich eine bestimmte Stichprobenverteilung bei kleinerem \(n\) ebenfalls annähernd normalverteilt³. Dann sind die Quantile bekannt und man kann die Streuung der Schätzers, \({\sigma }_{\bar {x}}\), z.B. für den Mittelwert, einfach schätzen:

\[{\displaystyle {\sigma }_{\bar {x}}\ ={\frac {\sigma }{\sqrt {n}}}}\]

7.5.3 Resampling passiert im Train-Sample

Wichtig zu beachten ist, dass die Resampling nur im Train-Sample stattfindet. Das Test-Sample bleibt unangerührt. Dieser Sachverhalt ist in Abb. 7.6, aus Silge and Kuhn (2022), illustriert.

Abbildung 7.6: Resampling im Train-, nicht im Test-Sample

Wie in Abb. 7.6 dargestellt, wird das Modell im Analyse-Sample berechnet (gefittet), und im Assessment-Sample auf Modellgüte hin überprüft.

Die letztliche Modellgüte ist dann die Zusammenfassung (Mittelwert) der einzelnen Resamples.

7.5.4 Andere Illustrationen

Es gibt eine Reihe vergleichbarer Illustrationen in anderen Büchern:

• Timbers, Campbell & Lee, 2022, Kap. 6
• Silge & Kuhn, 2022, Abb. 10.1
• Silge & Kuhn, 2022, Abb. 10.2
• Silge & Kuhn, 2022, Abb. 10.3
• James, Witten, hastie & Tishirani, 2021, Abb. 5.3

7.8.1 CV definieren

So kann man eine einfache v-fache Kreuzvalidierung in Tidymodels auszeichnen:

set.seed(2453)
ames_folds <- vfold_cv(ames_train, strata = "Sale_Price")
ames_folds

## #  10-fold cross-validation using stratification 
## # A tibble: 10 × 2
##    splits             id    
##    <list>             <chr> 
##  1 <split [1976/221]> Fold01
##  2 <split [1976/221]> Fold02
##  3 <split [1976/221]> Fold03
##  4 <split [1976/221]> Fold04
##  5 <split [1977/220]> Fold05
##  6 <split [1977/220]> Fold06
##  7 <split [1978/219]> Fold07
##  8 <split [1978/219]> Fold08
##  9 <split [1979/218]> Fold09
## 10 <split [1980/217]> Fold10

Werfen wir einen Blick in die Spalte splits, erste Zeile:

ames_folds %>% pluck(1, 1)

## <Analysis/Assess/Total>
## <1976/221/2197>

Möchte man die Defaults vpn vfold_cv wissen, schaut man in der Hilfe nach: ?vfold_cv:

vfold_cv(data, v = 10, repeats = 1, strata = NULL, breaks = 4, pool = 0.1, ...)

Probieren wir \(v=5\) und \(r=2\):

ames_folds_rep <- vfold_cv(ames_train, 
                           strata = "Sale_Price", 
                           v = 5,
                           repeats = 2)
ames_folds_rep

## #  5-fold cross-validation repeated 2 times using stratification 
## # A tibble: 10 × 3
##    splits             id      id2  
##    <list>             <chr>   <chr>
##  1 <split [1756/441]> Repeat1 Fold1
##  2 <split [1757/440]> Repeat1 Fold2
##  3 <split [1757/440]> Repeat1 Fold3
##  4 <split [1758/439]> Repeat1 Fold4
##  5 <split [1760/437]> Repeat1 Fold5
##  6 <split [1756/441]> Repeat2 Fold1
##  7 <split [1757/440]> Repeat2 Fold2
##  8 <split [1757/440]> Repeat2 Fold3
##  9 <split [1758/439]> Repeat2 Fold4
## 10 <split [1760/437]> Repeat2 Fold5

7.8.2 Resamples fitten

Hat unser Computer mehrere Rechenkerne, dann können wir diese nutzen und die Berechnungen beschleunigen. Im Standard wird sonst nur ein Kern verwendet.

mycores <- parallel::detectCores(logical = FALSE)
mycores

## [1] 4

Auf Unix/MacOC-Systemen kann man dann die Anzahl der parallen Kerne so einstellen:

library(doMC)
registerDoMC(cores = mycores)

So, und jetzt fitten wir die Resamples und trachten die Modellgüte in den Resamples:

ames_resamples_fit <- 
  ames_wflow %>% 
  fit_resamples(ames_folds)

 ames_resamples_fit %>%
  collect_metrics()

## # A tibble: 2 × 6
##   .metric .estimator   mean     n std_err .config             
##   <chr>   <chr>       <dbl> <int>   <dbl> <chr>               
## 1 rmse    standard   0.0928    10 0.00187 Preprocessor1_Model1
## 2 rsq     standard   0.722     10 0.00864 Preprocessor1_Model1

Natürlich interessiert uns primär die Modellgüte im Test-Sample:

final_ames <-
  last_fit(ames_wflow, data_split)

final_ames %>% 
  collect_metrics()

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard       0.103 Preprocessor1_Model1
## 2 rsq     standard       0.678 Preprocessor1_Model1

7.9.1 Tuning auszeichnen

In der Modellspezifikation des Modells können wir mit tune() auszeichnen, welche Parameter wir tunen möchten. Wir könenn

knn_model <-
  nearest_neighbor(
    mode = "regression",
    neighbors = tune()
  ) %>%
  set_engine("kknn")

Wir können dem Tuningparameter auch einen Namen (ID/Laben) geben, z.B. “K”:

knn_model <-
  nearest_neighbor(
    mode = "regression",
    neighbors = tune("K")
  ) %>%
  set_engine("kknn")

7.9.2 Grid Search vs. Iterative Search

Im K-Nächste-Nachbarn-Modell ist der vorhergesagt Wert, \(\hat{y}\) für eine neue Beobachtung \(x_0\) wie folgt definiert:

\[ \hat y = \frac{1}{K}\sum_{\ell = 1}^K x_\ell^*, \]

wobei \(K\) die Anzahl der zu berücksichtigen nächsten Nachbarn darstellt und \(x_\ell^*\) die Werte dieser berücksichtiggten Nachbarn.

Die Wahl von \(K\) hat einen gewaltigen Einfluss auf die Vorhersagen und damit auf die Vorhersagegüte. Allerdings wird \(K\) nicht vom Modell geschätzt. Es liegt an den Nutzi, diesen Wert zu wählen.

Parameter dieser Art (die von den Nutzi zu bestimmen sind, nicht vom Algorithmus), nennt man Tuningparameter.

Abbildung 7.9 aus Silge and Kuhn (2022) stellt exemplarisch dar, welchen großen Einfluss die Wahl des Werts eines Tuningparameters auf die Vorhersagen eines Modells haben.

Abbildung 7.9: Overfitting als Funktion der Modellparameter und insofern als Problem de Wahl der Tuningparameter

Aber wie wählt man “gute” Werte der Tuningparater? Zwei Ansätze, grob gesprochen, bieten sich an.

1. Grid Search: Probiere viele Werte aus und schaue, welcher der beste ist. Dabei musst du hoffen, dass du die Werte erwischt, die nicht nur im Train-, sondern auch im Test-Sample gut funktionieren werden.

2. Iterative Search: Wenn du einen Wert eines Tuningparameters hast, nutze diesen, um intelligenter einen neuen Wert eines Tuningparameters zu finden.

Der Unterschied beider Ansätze ist in Silge and Kuhn (2022) wie in Abb. 7.10 dargestellt.

Abbildung 7.10: Links: Grid Search. Rechts: Iterative Search2

In tidymodels kann man mit tune() angeben, dass man einen bestimmten Parameter tunen möchte. tidymodels führt das dann ohne weiteres Federlesens für uns durch.

7.10.0.1 Tuningparameter betrachten

Möchte man wissen, welche und wie viele Tuningparameter tidymodels in einem Modell berücksichtigt, kann man extract_parameter_set_dials() aufrufen:

extract_parameter_set_dials(knn_model)

## Collection of 1 parameters for tuning
## 
##  identifier      type    object
##           K neighbors nparam[+]

Die Ausgabe informiert uns, dass es nur einen Tuningparameter gibt in diesem Modell und dass der Name (Label, ID) des Tuningparameters “K” ist. Außerdem sollen die Anzahl der Nachbarn getunt werden. Der Tuningparameter ist numerisch; das sieht man an nparam[+].

Schauen wir uns mal an, auf welchen Wertebereich tidymodels den Parameter \(K\) begrenzt hat:

knn_model %>% 
  extract_parameter_dials("K")

## # Nearest Neighbors (quantitative)
## Range: [1, 15]

Aktualisieren wir mal unseren Workflow entsprechend:

ames_wflow <-
  ames_wflow %>% 
  update_model(knn_model)

Wir können auch Einfluss nehmen und angeben, dass die Grenzen des Wertebereichs zwischen 1 und 50 liegen soll (für den Tuningparameter neighbors):

ames_set <-
  extract_parameter_set_dials(ames_wflow) %>%
  update(K = neighbors(c(1, 50)))

ames_set

## Collection of 1 parameters for tuning
## 
##  identifier      type    object
##           K neighbors nparam[+]

7.10.1 Datenabhängige Tuningparameter

Manche Tuningparameter kann man nur bestimmen, wenn man den Datensatz kennt. So ist die Anzahl der Prädiktoren, mtry in einem Random-Forest-Modell sinnvollerweise als Funktion der Prädiktorenzahl zu wählen. Der Workflow kennt aber den Datensatz nicht. Daher muss der Workflow noch “finalisiert” oder “aktualisiert” werden, um den Wertebereich (Unter- und Obergrenze) eines Tuningparameters zu bestimmen.

Wenn wir im Rezept aber z.B. die Anzahl der Prädiktoren verändert haben, möchten wir die Grenzen des Wertebereichs für mtry (oder andere Tuningparameter) vielleicht nicht händisch, “hartverdrahtet” selber bestimmen, sondern lieber den Computer anweisen, und sinngemäß sagen: “Warte mal mit der Bestimmung der Werte der Tuningparameter, bis du den Datensatz bzw. dessen Dimensionen kennst. Merk dir, dass du, wenn du den Datensatz kennst, die Werte des Tuningparameter noch ändern musst. Und tu das dann auch.” Dazu später mehr.

ames_set <-
  workflow() %>% 
  add_model(knn_model) %>% 
  add_recipe(ames_rec) %>% 
  extract_parameter_set_dials() %>% 
  finalize(ames_train)

7.10.2 Modelle mit Tuning berechnen

Nachdem wir die Tuningwerte bestimmt haben, können wir jetzt das Modell berechnen: Für jeden Wert des Tuningparameters wird ein Modell berechnet:

ames_grid_search <-
  tune_grid(
    ames_wflow,
    resamples = ames_folds
  )
ames_grid_search

## # Tuning results
## # 10-fold cross-validation using stratification 
## # A tibble: 10 × 4
##    splits             id     .metrics          .notes          
##    <list>             <chr>  <list>            <list>          
##  1 <split [1976/221]> Fold01 <tibble [16 × 5]> <tibble [0 × 3]>
##  2 <split [1976/221]> Fold02 <tibble [16 × 5]> <tibble [0 × 3]>
##  3 <split [1976/221]> Fold03 <tibble [16 × 5]> <tibble [0 × 3]>
##  4 <split [1976/221]> Fold04 <tibble [16 × 5]> <tibble [0 × 3]>
##  5 <split [1977/220]> Fold05 <tibble [16 × 5]> <tibble [0 × 3]>
##  6 <split [1977/220]> Fold06 <tibble [16 × 5]> <tibble [0 × 3]>
##  7 <split [1978/219]> Fold07 <tibble [16 × 5]> <tibble [0 × 3]>
##  8 <split [1978/219]> Fold08 <tibble [16 × 5]> <tibble [0 × 3]>
##  9 <split [1979/218]> Fold09 <tibble [16 × 5]> <tibble [0 × 3]>
## 10 <split [1980/217]> Fold10 <tibble [16 × 5]> <tibble [0 × 3]>

Im Default berechnet tiymodels 10 Kandidatenmodelle.

Die Spalte .metrics beinhaltet die Modellgüte für jedes Kandidatenmodell.

ames_grid_search %>% 
  collect_metrics()

## # A tibble: 16 × 7
##        K .metric .estimator   mean     n std_err .config             
##    <int> <chr>   <chr>       <dbl> <int>   <dbl> <chr>               
##  1     2 rmse    standard   0.103     10 0.00213 Preprocessor1_Model1
##  2     2 rsq     standard   0.662     10 0.0112  Preprocessor1_Model1
##  3     4 rmse    standard   0.0950    10 0.00188 Preprocessor1_Model2
##  4     4 rsq     standard   0.708     10 0.00916 Preprocessor1_Model2
##  5     6 rmse    standard   0.0912    10 0.00189 Preprocessor1_Model3
##  6     6 rsq     standard   0.732     10 0.00842 Preprocessor1_Model3
##  7     7 rmse    standard   0.0900    10 0.00192 Preprocessor1_Model4
##  8     7 rsq     standard   0.740     10 0.00829 Preprocessor1_Model4
##  9     9 rmse    standard   0.0883    10 0.00201 Preprocessor1_Model5
## 10     9 rsq     standard   0.752     10 0.00827 Preprocessor1_Model5
## 11    11 rmse    standard   0.0872    10 0.00211 Preprocessor1_Model6
## 12    11 rsq     standard   0.761     10 0.00845 Preprocessor1_Model6
## 13    13 rmse    standard   0.0865    10 0.00217 Preprocessor1_Model7
## 14    13 rsq     standard   0.767     10 0.00848 Preprocessor1_Model7
## 15    15 rmse    standard   0.0861    10 0.00221 Preprocessor1_Model8
## 16    15 rsq     standard   0.772     10 0.00850 Preprocessor1_Model8

Das können wir uns einfach visualisieren lassen:

autoplot(ames_grid_search)

Auf Basis dieser Ergebnisse könnte es Sinn machen, noch größere Werte für \(K\) zu überprüfen.

7.10.3 Vorhersage im Test-Sample

Welches Modellkandidat war jetzt am besten?

show_best(ames_grid_search)

## # A tibble: 5 × 7
##       K .metric .estimator   mean     n std_err .config             
##   <int> <chr>   <chr>       <dbl> <int>   <dbl> <chr>               
## 1    15 rmse    standard   0.0861    10 0.00221 Preprocessor1_Model8
## 2    13 rmse    standard   0.0865    10 0.00217 Preprocessor1_Model7
## 3    11 rmse    standard   0.0872    10 0.00211 Preprocessor1_Model6
## 4     9 rmse    standard   0.0883    10 0.00201 Preprocessor1_Model5
## 5     7 rmse    standard   0.0900    10 0.00192 Preprocessor1_Model4

Wählen wir jetzt mal das beste Modell aus (im Sinne des Optimierungskriteriusms):

select_best(ames_grid_search)

## # A tibble: 1 × 2
##       K .config             
##   <int> <chr>               
## 1    15 Preprocessor1_Model8

Ok, notieren wir uns die Kombination der Tuningparameterwerte im besten Kandiatenmodell. In diesem Fall hat das Modull nur einen Tuningparameter:

ames_knn_best_params <-
  tibble(K = 15)

Unser Workflow weiß noch nicht, welche Tuningparameterwerte am besten sind:

ames_wflow

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: nearest_neighbor()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 4 Recipe Steps
## 
## • step_log()
## • step_other()
## • step_dummy()
## • step_zv()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## K-Nearest Neighbor Model Specification (regression)
## 
## Main Arguments:
##   neighbors = tune("K")
## 
## Computational engine: kknn

neighbors = tune("K") sagt uns, dass er diesen Parameter tunen will. Das haben wir jetzt ja erledigt. Wir wollen für das Test-Sample nur noch einen Wert, eben aus dem besten Kandidatenmodell, verwenden:

ames_final_wflow <-
  ames_wflow %>% 
  finalize_workflow(ames_knn_best_params)

ames_final_wflow

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: nearest_neighbor()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 4 Recipe Steps
## 
## • step_log()
## • step_other()
## • step_dummy()
## • step_zv()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## K-Nearest Neighbor Model Specification (regression)
## 
## Main Arguments:
##   neighbors = 15
## 
## Computational engine: kknn

Wie man sieht, steht im Workflow nichts mehr von Tuningparameter.

Wir können jetzt das ganze Train-Sample fitten, also das Modell auf das ganze Train-Sample anwenden - nicht nur auf ein Analysis-Sample. Und mit den dann resultierenden Modellkoeffizienten sagen wir das TestSample vorher:

final_ames_knn_fit <-
  last_fit(ames_final_wflow, data_split)

final_ames_knn_fit

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits             id               .metrics .notes   .predictions .workflow 
##   <list>             <chr>            <list>   <list>   <list>       <list>    
## 1 <split [2197/733]> train/test split <tibble> <tibble> <tibble>     <workflow>

Holen wir uns die Modellgüte:

collect_metrics(final_ames_knn_fit)

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard      0.0951 Preprocessor1_Model1
## 2 rsq     standard      0.736  Preprocessor1_Model1

8.1.1 Vorbereitung

• Frischen Sie Ihr Wissen zur logistischen Regression auf bzw. machen Sie sich mit den Grundlagen des Verfahrens vertraut.

8.1.2 Lernziele

• Sie verstehen den Zusammenhang von linearen und logistischen Modellen
• Sie können die logistische Regression mit Methoden von tidymodels anwenden

8.1.3 Literatur

• Rhys, Kap. 4

8.4.1 Lineare Modelle running wild

Wie groß ist die Wahrscheinlichkeit für eine manuelle Schaltung …

• … bei mpg_z = -2?

predict(m81, newdata = data.frame(iv = -2))

##          1 
## -0.1923719

\(Pr(\hat{y})<0\) macht keinen Sinn. ⚡

• … bei mpg_z = +2?

predict(m81, newdata = data.frame(iv = +2))

##        1 
## 1.004872

\(Pr(\hat{y})>1\) macht keinen Sinn. ⚡

Schauen Sie sich mal die Vorhersage an für mpg_z=5 🤯

8.4.2 Wir müssen die Regressionsgerade umbiegen

… wenn der vorhergesagte Wert eine Wahrscheinlichkeit, \(p_i\), ist.

Die schwarze Gerade verlässt den Wertebereich der Wahrscheinlichkeit. Die blaue Kurve, \(\mathcal{f}\), bleibt im erlaubten Bereich, \(Pr(y) \in [0,1]\). Wir müssen also die linke oder die rechte Seite des linearen Modells transformieren: \(p_i = f(\alpha + \beta \cdot x)\) bzw.:

\(f(p) = \alpha + \beta \cdot x\)

\(\mathcal{f}\) nennt man eine Link-Funktion.

8.4.3 Verallgemeinerte lineare Modelle zur Rettung

Für metrische AV mit theoretisch unendlichen Grenzen des Wertebereichs haben wir bisher eine Normalverteilung verwendet:

\[y_i \sim \mathcal{N}(\mu_i, \sigma)\]

Dann ist die Normalverteilung eine voraussetzungsarme Wahl (maximiert die Entropie).

Aber wenn die AV binär ist bzw. Häufigkeiten modelliert, braucht man eine Variable die nur positive Werte zulässt.

Diese Verallgemeinerung des linearen Modells bezeichnet man als verallgemeinertes lineares Modell (generalized linear model, GLM).

Im Falle einer binären (bzw. dichotomen) AV liegt eine bestimmte Form des GLM vor, die man als logistische Regression bezeichnet.

8.6.1 tidymodels, m83

Berechnen wir jetzt ein lineares Modell für die AV Survived mit dem Geschlecht als Pädiktor:

d <-
  titanic_train %>% 
  filter(Fare > 0) %>% 
  mutate(iv = log(Fare),
         dv = factor(Survived))

Die Faktorstufen, genannt levels von Survived sind:

levels(d$dv)

## [1] "0" "1"

Und zwar genau in dieser Reihenfolge.

8.8.1 Wahrscheinlichkeit in Odds

Probieren wir Folgendes: Rechnen wir die Wahrscheinlichkeit zu Überlegen für \(y\), kurz \(p\), in Odds (Chancen) um.

\(odds = \frac{p}{1-p}\)

In R:

odds <- 0.38 / 0.62
odds

## [1] 0.6129032

Bildlich gesprochen sagen die Odds: für 38 Menschen, die überlebt haben, kommen (ca.) 62 Menschen, die nicht überlebt haben, s. Abb. 8.3.

Abbildung 8.3: Odds: 38 zu 62

Plotten wir die Odds als Funktion der UV, s. Abb. 8.4.

Abbildung 8.4: Odds als Funktion der UV

Wir sind noch nicht am Ziel; die Variable ist noch nicht “richtig gebogen”.

8.8.2 Von Odds zu Log-Odds

Wenn wir jetzt den Logarithmus berechnen der Log-Odds bekommen wir eine “brav gebogenen” Funktion, die Log-Odds, \(\mathcal{L}\), als Funktion der UV, s. Abb. 8.5.

\[\mathcal{L} = log (odds) = log \left(\frac{p}{1-p}\right)\]

Abbildung 8.5: Logit als Funktion der UV

Linear!

Es gilt also:

\[\text {log-odds} = b_0 + b_1x\]

Log-Odds (Log-Odds) bezeichnet man auch als Logits.

8.10.1 Grenzwert wechseln

Im Standard wird 50% als Grenzwert für die vorhergesagte Klasse \(c\) genommen:

• wenn \(p <= .5 \rightarrow c = 0\)
• wenn \(p > .5 \rightarrow c = 1\)

Man kann aber den Grenzwert beliebig wählen, um Kosten-Nutzen-Abwägungen zu optimieren; mehr dazu findet sich z.B. hier.

8.11.1 Logit

\((0,1) \rightarrow (-\infty, +\infty)\)

Praktisch, um Wahrscheinlichkeit zu modellieren.

\[p \rightarrow \fbox{logit} \rightarrow \alpha + \beta x\]

8.11.2 Inv-Logit

\((-\infty, +\infty) \rightarrow (0,1)\)

Praktisch, um in Wahrscheinlichkeiten umzurechnen.

\[p \leftarrow \fbox{inv-logit} \leftarrow \alpha + \beta x\]

8.12.1 Die Koeffizienten sind schwer zu interpretieren

• In der logistischen Regression gilt nicht mehr, dass eine konstante Veränderung in der UV mit einer konstanten Veränderung in der AV einhergeht.
• Stattdessen geht eine konstante Veränderung in der UV mit einer konstanten Veränderung im Logit der AV einher.
• Beim logistischen Modell hier gilt, dass in der Nähe von \(x=0\) die größte Veränderung in \(p\) von statten geht; je weiter weg von \(x=0\), desto geringer ist die Veränderung in \(p\).

8.12.2 Logits vs. Wahrscheinlichkeiten

konvert_logits <-
  tibble(
    logit = c( -10, -3, 
              -2, -1, -0.5, -.25, 
              0, 
              .25, .5, 1, 2, 
              3, 10),
    p = rstanarm::invlogit(logit)
  )  %>% 
  gt() %>% 
  fmt_number(everything(), decimals = 2)

9.11.1 Initiale Datenaufteilung

library(tidymodels)

data("mtcars")

set.seed(42)  # Reproduzierbarkeit
d_split <- initial_split(mtcars, strata = mpg)

## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 1 breaks instead.

## Warning: Too little data to stratify.
## • Resampling will be unstratified.

d_train <- training(d_split)
d_test <- testing(d_split)

Die Warnung zeigt uns, dass der Datensatz sehr klein ist; stimmt. Ignorieren wir hier einfach.

Wie man auf der Hilfeseite der Funktion sieht, wird per Voreinstellung 3/1 aufgeteilt, also 75% in das Train-Sample, 25% der Daten ins Test-Sample.

Bei \(n=32\) finden also 8 Autos ihren Weg ins Test-Sample und die übrigen 24 ins Train-Sample. Bei der kleinen Zahl könnte man sich (berechtigterweise) fragen, ob es Sinn macht, die spärlichen Daten noch mit einem Test-Sample weiter zu dezimieren. Der Einwand ist nicht unberechtigt, allerdings zieht der Verzicht auf ein Test-Sample andere Probleme, Overfitting namentlich, nach sich.

9.11.2 Kreuzvalidierung definieren

d_cv <- vfold_cv(d_train, strata = mpg, repeats = 5, v = 5) 
d_cv

## #  5-fold cross-validation repeated 5 times using stratification 
## # A tibble: 25 × 3
##    splits         id      id2  
##    <list>         <chr>   <chr>
##  1 <split [19/5]> Repeat1 Fold1
##  2 <split [19/5]> Repeat1 Fold2
##  3 <split [19/5]> Repeat1 Fold3
##  4 <split [19/5]> Repeat1 Fold4
##  5 <split [20/4]> Repeat1 Fold5
##  6 <split [19/5]> Repeat2 Fold1
##  7 <split [19/5]> Repeat2 Fold2
##  8 <split [19/5]> Repeat2 Fold3
##  9 <split [19/5]> Repeat2 Fold4
## 10 <split [20/4]> Repeat2 Fold5
## # … with 15 more rows

Die Defaults (Voreinstellungen) der Funktion vfold_cv() können, wie immer, auf der Hilfeseite der Funktion nachgelesen werden.

Da die Stichprobe sehr klein ist, bietet es sich an, eine kleine Zahl an Faltungen (folds) zu wählen. Bei 10 Faltungen beinhaltete eine Stichprobe gerade 10% der Fälle in Train-Sample, also etwa … 2!

Zur Erinnerung: Je größer die Anzahl der Repeats, desto genauer schätzen wir die Modellgüte.

9.11.3 Rezept definieren

Hier ein einfaches Rezept:

recipe1 <-
  recipe(mpg ~ ., data = d_train) %>% 
  step_impute_knn() %>% 
  step_normalize() %>% 
  step_dummy() %>% 
  step_other(threshold = .1)

9.11.4 Modell definieren

tree_model <-
  decision_tree(
    cost_complexity = tune(),
    tree_depth = tune(),
    min_n = tune()
  ) %>% 
  set_engine("rpart") %>% 
  set_mode("regression")

Wenn Sie sich fragen, woher Sie die Optionen für die Tuningparameter wissen sollen: Schauen Sie mal in die Hilfeseite des Pakets {{dials}}; das Paket ist Teil von Tidymodels.

Die Berechnung des Modells läuft über das Paket {{rpart}}, was wir durch set_engine() festgelegt haben.

Der Parameter Cost Complexity, \(C_p\) oder manchmal auch mit \(\alpha\) bezeichnet, hat einen typischen Wertebereich von \(10^{-10}\) bis \(10^{-1}\):

cost_complexity()

## Cost-Complexity Parameter (quantitative)
## Transformer: log-10 [1e-100, Inf]
## Range (transformed scale): [-10, -1]

Hier ist der Wert in Log-Einheiten angegeben. Wenn Sie sich fragen, woher Sie das bitteschön wissen sollen: Naja, es steht auf der Hilfeseite 😄.

Unser Modell ist also so definiert:

tree_model

## Decision Tree Model Specification (regression)
## 
## Main Arguments:
##   cost_complexity = tune()
##   tree_depth = tune()
##   min_n = tune()
## 
## Computational engine: rpart

Mit tune() weist man den betreffenden Parameter als “zu tunen” aus - gute Werte sollen durch Ausprobieren während des Berechnens bestimmt werden. Genauer gesagt soll das Modell für jeden Wert (oder jede Kombination an Werten von Tuningparametern) berechnet werden.

Eine Kombination an Tuningparameter-Werten, die ein Modell spezifizieren, sozusagen erst “fertig definieren”, nennen wir einen Modellkandidaten.

Definieren wir also eine Tabelle (grid) mit Werten, die ausprobiert, “getuned” werden sollen. Wir haben oben dre Tuningparameter bestimmt. Sagen wir, wir hätten gerne jeweils 5 Werte pro Parameter.

tree_grid <-
  grid_regular(
    cost_complexity(),
    tree_depth(),
    min_n(),
    levels = 4
  )

Für jeden Parameter sind Wertebereiche definiert; dieser Wertebereich wird gleichmäßig (daher grid regular) aufgeteilt; die Anzahl der verschiedenen Werte pro Parameter wird druch levels gegeben.

Mehr dazu findet sich auf der Hilfeseite zu grid_regular().

Wenn man die alle miteinander durchprobiert, entstehen \(4^3\) Kombinationen, also Modellkandidaten.

Allgemeiner gesagt sind das bei \(n\) Tuningparametern mit jeweils \(m\) verschiedenen Werten \(m^n\) Möglichkeiten, spricht Modellkandidaten. Um diesen Faktor erhöht sich die Rechenzeit im Vergleich zu einem Modell ohne Tuning. Man sieht gleich, dass die Rechenzeit schnell unangenehm lang werden kann.

Entsprechend hat unsere Tabelle diese Zahl an Zeilen. Jede Zeile definiert einen Modellkandidaten, also eine Berechnung des Modells.

dim(tree_grid)

## [1] 64  3

head(tree_grid)

## # A tibble: 6 × 3
##   cost_complexity tree_depth min_n
##             <dbl>      <int> <int>
## 1    0.0000000001          1     2
## 2    0.0000001             1     2
## 3    0.0001                1     2
## 4    0.1                   1     2
## 5    0.0000000001          5     2
## 6    0.0000001             5     2

Man beachte, dass außer Definitionen bisher nichts passiert ist – vor allem haben wir noch nichts berechnet. Sie scharren mit den Hufen? Wollen endlich loslegen? Also gut.

9.11.5 Workflow definieren

Fast vergessen: Wir brauchen noch einen Workflow.

tree_wf <-
  workflow() %>% 
  add_model(tree_model) %>% 
  add_recipe(recipe1)

9.11.6 Modell tunen und berechnen

Achtung: Das Modell zu berechnen kann etwas dauern. Es kann daher Sinn machen, das Modell abzuspeichern, so dass Sie beim erneuten Durchlaufen nicht nochmal berechnen müssen, sondern einfach von der Festplatte laden können; das setzt natürlich voraus, dass sich am Modell nichts geändert hat.

doParallel::registerDoParallel()  # mehrere Kerne parallel nutzen

set.seed(42)
tic()  # Stoppuhr an
trees_tuned <-
  tune_grid(
    object = tree_wf,
    grid = tree_grid,
    resamples = d_cv
  )
toc()  # Stoppuhr aus

Es bietet sich in dem Fall an, das Ergebnis-Objekt als R Data serialized (rds) abzuspeichern:

write_rds(trees_tuned, "objects/trees1.rds")

Bzw. so wieder aus der RDS-Datei zu importieren:

trees_tuned <- read_rds("objects/trees1.rds")

Hier oder hier kann man einiges zum Unterschied einer RDS-Datei vs. einer “normalen” R-Data-Datei nachlesen. Wenn man möchte 😉.

trees_tuned

## # Tuning results
## # 5-fold cross-validation repeated 5 times using stratification 
## # A tibble: 25 × 5
##    splits         id      id2   .metrics           .notes           
##    <list>         <chr>   <chr> <list>             <list>           
##  1 <split [19/5]> Repeat1 Fold1 <tibble [128 × 7]> <tibble [32 × 3]>
##  2 <split [19/5]> Repeat1 Fold2 <tibble [128 × 7]> <tibble [32 × 3]>
##  3 <split [19/5]> Repeat1 Fold3 <tibble [128 × 7]> <tibble [32 × 3]>
##  4 <split [19/5]> Repeat1 Fold4 <tibble [128 × 7]> <tibble [32 × 3]>
##  5 <split [20/4]> Repeat1 Fold5 <tibble [128 × 7]> <tibble [32 × 3]>
##  6 <split [19/5]> Repeat2 Fold1 <tibble [128 × 7]> <tibble [32 × 3]>
##  7 <split [19/5]> Repeat2 Fold2 <tibble [128 × 7]> <tibble [32 × 3]>
##  8 <split [19/5]> Repeat2 Fold3 <tibble [128 × 7]> <tibble [32 × 3]>
##  9 <split [19/5]> Repeat2 Fold4 <tibble [128 × 7]> <tibble [32 × 3]>
## 10 <split [20/4]> Repeat2 Fold5 <tibble [128 × 7]> <tibble [32 × 3]>
## # … with 15 more rows
## 
## There were issues with some computations:
## 
##   - Warning(s) x320: 27 samples were requested but there were 19 rows in the data. 19 ...   - Warning(s) x80: 27 samples were requested but there were 19 rows in the data. 19 ...   - Warning(s) x320: 27 samples were requested but there were 19 rows in the data. 19 ...   - Warning(s) x80: 27 samples were requested but there were 19 rows in the data. 19 ...   - Warning(s) x4: 27 samples were requested but there were 19 rows in the data. 19 ...
## 
## Use `collect_notes(object)` for more information.

Die Warnhinweise kann man sich so ausgeben lassen:

collect_notes(trees_tuned)

## # A tibble: 804 × 5
##    id      id2   location                      type    note                     
##    <chr>   <chr> <chr>                         <chr>   <chr>                    
##  1 Repeat1 Fold1 preprocessor 1/1, model 33/64 warning 27 samples were requeste…
##  2 Repeat1 Fold1 preprocessor 1/1, model 34/64 warning 27 samples were requeste…
##  3 Repeat1 Fold1 preprocessor 1/1, model 35/64 warning 27 samples were requeste…
##  4 Repeat1 Fold1 preprocessor 1/1, model 36/64 warning 27 samples were requeste…
##  5 Repeat1 Fold1 preprocessor 1/1, model 37/64 warning 27 samples were requeste…
##  6 Repeat1 Fold1 preprocessor 1/1, model 38/64 warning 27 samples were requeste…
##  7 Repeat1 Fold1 preprocessor 1/1, model 39/64 warning 27 samples were requeste…
##  8 Repeat1 Fold1 preprocessor 1/1, model 40/64 warning 27 samples were requeste…
##  9 Repeat1 Fold1 preprocessor 1/1, model 41/64 warning 27 samples were requeste…
## 10 Repeat1 Fold1 preprocessor 1/1, model 42/64 warning 27 samples were requeste…
## # … with 794 more rows

Wie gesagt, in diesem Fall war die Stichprobengröße sehr klein.

9.11.7 Modellgüte evaluieren

collect_metrics(trees_tuned)

## # A tibble: 128 × 9
##    cost_complexity tree_depth min_n .metric .estimator  mean     n std_err
##              <dbl>      <int> <int> <chr>   <chr>      <dbl> <int>   <dbl>
##  1    0.0000000001          1     2 rmse    standard   3.46     25  0.223 
##  2    0.0000000001          1     2 rsq     standard   0.666    21  0.0385
##  3    0.0000001             1     2 rmse    standard   3.46     25  0.223 
##  4    0.0000001             1     2 rsq     standard   0.666    21  0.0385
##  5    0.0001                1     2 rmse    standard   3.46     25  0.223 
##  6    0.0001                1     2 rsq     standard   0.666    21  0.0385
##  7    0.1                   1     2 rmse    standard   3.46     25  0.223 
##  8    0.1                   1     2 rsq     standard   0.666    21  0.0385
##  9    0.0000000001          5     2 rmse    standard   2.62     25  0.265 
## 10    0.0000000001          5     2 rsq     standard   0.863    25  0.0279
## # … with 118 more rows, and 1 more variable: .config <chr>

Praktischerweise gibt es eine Autoplot-Funktion, um die besten Modellparameter auszulesen:

autoplot(trees_tuned)

9.11.8 Bestes Modell auswählen

Aus allen Modellkandidaten wählen wir jetzt das beste Modell aus:

select_best(trees_tuned)

## # A tibble: 1 × 4
##   cost_complexity tree_depth min_n .config              
##             <dbl>      <int> <int> <chr>                
## 1          0.0001          5     2 Preprocessor1_Model07

Mit diesem besten Kandidaten definieren wir jetzt das “finale” Modell, wir “finalisieren” das Modell mit den besten Modellparametern:

tree_final <-
  finalize_model(tree_model, parameters = select_best(trees_tuned))

tree_final

## Decision Tree Model Specification (regression)
## 
## Main Arguments:
##   cost_complexity = 1e-04
##   tree_depth = 5
##   min_n = 2
## 
## Computational engine: rpart

Hier ist, unser finaler Baum 🌳.

Schließlich updaten wir mit dem finalen Baum noch den Workflow:

final_wf <-
  tree_wf %>% 
  update_model(tree_final)

9.11.9 Final Fit

Jetzt fitten wir dieses Modell auf das ganze Train-Sample und predicten auf das Test-Sample:

tree_fit_final <-
  final_wf %>% 
  last_fit(d_split)

tree_fit_final

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits         id               .metrics .notes   .predictions     .workflow 
##   <list>         <chr>            <list>   <list>   <list>           <list>    
## 1 <split [24/8]> train/test split <tibble> <tibble> <tibble [8 × 4]> <workflow>

collect_metrics(tree_fit_final)

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard       3.93  Preprocessor1_Model1
## 2 rsq     standard       0.683 Preprocessor1_Model1

Voilà: Die Modellgüte für das Test-Sample: Im Schnitt liegen wir ca. 4 Meilen daneben mit unseren Vorhersagen, wenn wir RMSE mal so locker interpretieren wollen.

In der Regel ist übrigens RMSE interessanter als R-Quadrat, da R-Quadrat die Güte eines Korrelationsmusters vorhersagt, aber RMSE die Präzision der Vorhersage, also sozusagen die Kürze der Fehlerbalken.

9.11.10 Nur zum Spaß: Vergleich mit linearem Modell

Ein einfaches lineares Modell, was hätte das jetzt wohl für eine Modellgüte?

lm_model <-
  linear_reg()

lm_wf <-
  workflow() %>% 
  add_model(lm_model) %>% 
  add_recipe(recipe1)

tic()
lm_fit <-
  fit_resamples(
    lm_wf,
    resamples = d_cv
  )
toc()

## 7.662 sec elapsed

collect_metrics(lm_fit)

## # A tibble: 2 × 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 rmse    standard   4.16     25  0.362  Preprocessor1_Model1
## 2 rsq     standard   0.624    25  0.0587 Preprocessor1_Model1

lm_fit_final <- 
  last_fit(lm_wf, d_split)

Wie präzise ist die Vorhersage im Test-Sample?

collect_metrics(lm_fit_final)

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard       5.24  Preprocessor1_Model1
## 2 rsq     standard       0.434 Preprocessor1_Model1

Das lineare Modell schneidet etwas (deutlich?) schlechter ab als das einfache Baummodell.

Man beachte, dass die Modellgüte im Train-Samle höher ist als im Test-Sample (Overfitting).

10.1.1 Lernziele

• Sie können Algorithmen für Ensemble-Lernen erklären, d.i. Bagging, AdaBoost, XGBoost, Random Forest
• Sie wissen, anhand welche Tuningparamter man Overfitting bei diesen Algorithmen begrenzen kann
• Sie können diese Verfahren in R berechnen

10.1.2 Literatur

• Rhys, Kap. 8

10.1.3 Hinweise

• Nutzen Sie StackOverflow als Forum für Ihre Fragen - Hier ein Beispiel zu einer Fehlermeldung, die mir Kopfzerbrechen bereitete

10.6.1 Bootstrapping

Das erste baumbasierte Modell, was vorgestellt werden soll, basiert auf sog. Bootstrapping, ein Standardverfahren in der Statistik (James et al. 2021).

Bootstrapping ist eine Nachahmung für folgende Idee: Hätte man viele Stichproben aus der relevanten Verteilung, so könnte man z.B. die Genauigkeit eines Modells \(\hat{f}_{\bar{X}}\) zur Schätzung des Erwartungswertes \(\mu\) einfach dadurch bestimmen, indem man se berechnet, also die Streuung der Mitterwerte \(\bar{X}\) berechnet. Außerdem gilt, dass die Präzision der Schätzung des Erwartungswerts steigt mit steigendem Stichprobenumfang \(n\). Wir könnten also für jede der \(B\) Stichproben, \(b=1,\ldots, B\), ein (Baum-)Modell berechnen, \(\hat{f}^b\), und dann deren Vorhersagen aggregieren (zum Mittelwert oder Modalwert). Das kann man formal so darstellen (James et al. 2021):

\[\hat{f}_{\bar{X}} = \frac{1}{B}\sum_{b=1}^{B}\hat{f}^b\]

Mit diesem Vorgehen kann die Varianz des Modells \(\hat{f}_{\bar{X}}\) verringert werden; die Vorhersagegenauigkeit steigt.

Leider haben wir in der Regel nicht viele (\(B\)) Datensätze.

Daher “bauen” wir uns aus dem einzelnen Datensatz, der uns zur Verfügung steht, viele Datensätze. Das hört sich nach “too good to be true” an⁶ Weil es sich unglaubwürdig anhört, nennt man das entsprechende Verfahren (gleich kommt es!) auch “Münchhausen-Methode”, nach dem berühmten Lübgenbaron. Die Amerikaner ziehen sich übrigens nicht am Schopf aus dem Sumpf, sondern mit den Stiefelschlaufen (die Cowboys wieder), daher spricht man im Amerikanischen auch von der “Boostrapping-Methode”.

Diese “Pseudo-Stichproben” oder “Bootstrapping-Stichproben” sind aber recht einfach zu gewinnen.. Gegeben sei Stichprobe der Größe \(n\):

1. Ziehe mit Zurücklegen (ZmZ) aus der Stichprobe \(n\) Beobachtungen
2. Fertig ist die Bootstrapping-Stichprobe.

Abb. 10.1 verdeutlicht das Prinzip des ZMZ, d.h. des Bootstrappings. Wie man sieht, sind die Bootstrap-Stichproben (rechts) vom gleichen Umfang \(n\) wie die Originalstichprobe (links). Allerdins kommen nicht alle Fälle (in der Regel) in den “Boostrap-Beutel” (in bag), sondern einige Fälle werden oft mehrfach gezogen, so dass einige Fälle nicht gezogen werden (out of bag).

Abbildung 10.1: Bootstrapping: Der Topf links symbolisiert die Original-Stichprobe, aus der wir hier mehrere ZMZ-Stichproben ziehen (Rechts), dargestellt mit ‘in bag’

Man kann zeigen, dass ca. 2/3 der Fälle gezogen werden, bzw. ca. 1/3 nicht gezogen werden. Die nicht gezogenen Fälle nennt man auch out of bag (OOB).

Für die Entwicklung des Bootstrapping wurde der Autor, Bradley Efron, im Jahr 2018 mit dem internationalen Preis für Statistik ausgezeichnet;

“While statistics offers no magic pill for quantitative scientific investigations, the bootstrap is the best statistical pain reliever ever produced,” says Xiao-Li Meng, Whipple V. N. Jones Professor of Statistics at Harvard University.“

10.7.1 Variablenrelevanz

Man kann die Relevanz der Prädiktoren in einem Bagging-Modell auf mehrere Arten schätzen. Ein Weg (bei numerischer Prädiktion) ist, dass man die RSS-Verringerung, die durch Aufteilung anhand eines Prädiktors erzeugt wird, mittelt über alle beteiligten Bäume (Modelle). Bei Klassifikation kann man die analog die Reduktion des Gini-Wertes über alle Bäume mitteln und als Schätzwert für die Relevanz des Prädiktors heranziehen.

10.7.2 Out of Bag Vorhersagen

Da nicht alle Fälle der Stichprobe in das Modell einfließen (sondern nur ca. 2/3), kann der Rest der Fälle zur Vorhersage genutzt werden. Bagging erzeugt sozusagen innerhalb der Stichprobe selbständig ein Train- und ein Test-Sample. Man spricht von Out-of-Bag-Schätzung (OOB-Schätzung). Der OOB-Fehler (z.B. MSE bei numerischen Modellen und Genauigkeit bei nominalen) ist eine valide Schätzung des typischen Test-Sample-Fehlers.

Hat man aber Tuningparameter, so wird man dennoch auf die typische Train-Test-Aufteilung zurückgreifen, um Overfitting durch das Ausprobieren der Tuning-Kandidaten zu vermeiden (was sonst zu Zufallstreffern führen würde bei genügend vielen Modellkandidaten).